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Abstract 

A continuum theory to describe periodic ripple formation in planar graphene/boron nitride 
superlattices is formulated. Due to the lattice mismatch between the two materials, it is shown 
that flat superlattices are unstable with respect to ripple formation of appropriate wavelengths. 
A competition between bending energy and transverse stretching energy gives rise to an optimal 
ripple wavelength that depends on the superlattice pitch. The optimal wavelengths predicted by 
the continuum theory are in good agreement with atomic scale total energy calculations previously 
reported in Nandwana and Ertekin, Nano Lett. 15 1468 (2015). 

PACS numbers: 62.25.-g, 68.65.Pq, 68.60.Bs 
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I. INTRODUCTION 


The rippling, crumpling, and wrinkling of elastic materials are observed at many length 
scales, for example: crumpling of sheets of paper^’^, wrinkling of skin^, and the buckling 
of the earth’s crust due to viscous stresses from the underlying mantle^. The combination 
of geometry and physical models that capture the wavelength, amplitude, and curvature of 
large deformations often gives rise to surprising behavior, such as the famous problem of the 
wrinkling of an elastic sheet stretched in tension but clamped from contracting laterally at its 
ends^. When it comes to the deformation of naturally thin sheets, predicting the wavelength 
and amplitudes of key features has been a long-standing challenge®’^. The large deformations 
are governed by the Foppl-von Karman® equations. These non-linear equations reflect a 
subtle balance between stretching and bending energies in deformed membranes. With 
few exceptions, they are impossible to solve analytically, necessitating instead a variational 
approach and/or scaling analysis. 

Meanwhile, the recent development of two-dimensional and quasi- two-dimensional ma¬ 
terials such as graphene®, boron nitride^®, and transition metal dichalcogenides^^ now pro¬ 
vides an opportunity to observe the effects of elastic deformation at limiting length scales 
even down to a single atomic layer thick. For example, ripples are intrinsic to graphene 
sheets^^’^^, and are expected to strongly influence the electronic^^’^^, mechanicaF®’^^, and 
thermaF® properties. Control over the wavelength, amplitude, and orientation of one and 
two dimensional periodic ripples has been demonstrated^®, and it has been suggested that 
rippling can give rise to a negative thermal expansion coefficient^®. Also, thermal strain en¬ 
gineering and scanning tunneling microscopy of suspended nanographene membranes have 
suggested the breakdown of continuum theories when ripple wavelengths are close to the 
lattice constant^^. 

Our recent work focused on the prediction of the energetically favorable configurations of 
graphene/boron nitride planar superlattices^^. In these superlattices, a single atomic plane 
is comprised of perfectly ordered, alternating regions of boron nitride and graphene. As 
shown in Fig. la, due to the lattice mismatch between the two materials, the superlattice 
when flat is comprised of alternating regions of tension and compression along the direction 
parallel to the interface. While strain-relieving misfit dislocations are expected to relieve the 
lattice mismatch in super lattices that are constrained to remain flat, out-of-plane buckling 
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and ripple formation was found to also be a very effective strain relief avenue^^. 

In this work, we present a detailed continuum formulation and scaling analysis of rip¬ 
ple formation in these super lattices. Depending on the superlattice pitch H, the lattice 
mismatch renders the planar flat superlattices unstable with respect to the formation of a 
ripple network as shown in Fig. lb. The ripples exhibit larger amplitude in boron nitride to 
accommodate its larger length, but decay to smaller amplitude in graphene. We find that 
an optimal set of ripple wavelengths emerges from a competition between the ripple bending 
energy (favoring short wavelengths) and the transverse stretching energy due to ripple decay 
(favoring long wavelengths). Our continuum and scaling analysis is in excellent agreement 
with atomistic models of ripple formation^^, demonstrating that when properly formulated, 
continuum theories can hold quite nicely even at these length scales. Although the analysis 
presented here is specific to graphene/boron nitride superlattices, the framework is general 
and can be applied to other planar superlattices. 


II. PRELIMINARY DISCUSSION OF RIPPLING DISTORTIONS 

A. Total elastic energy of a deformed membrane. 

In general, the total elastic energy of a rippled sheet U = Us+Ub contains both stretching 
and bending contributions®: 


Us = J usdA = j 

f 2^^ijkFij^kldA , 

( 1 ) 

Ub — J UsdA = 

^ ^K(V‘^wfdA 

( 2 ) 


In the stretching texmUs^ Cijki are elastic constants and Cjj = ef^+Fij are total strains, where 
— \{diUj + djUi) is the distortion strain and Fij = \{diUkdjUk + diwdjw) is the bending 
strain. The indices i,j,k vary over x^y; dA = dxdy; and repeated indices are summed. 
The displacement field consists of in-plane u{x,y) = {ux{x,y),Uy{x,y)) and out-of-plane 
w{x,y) contributions. With large deflections, the nonlinear term Fij is non-negligible. The 
bending term Ub reflects contributions to the energy associated with changes to the gaussian 
curvature H ~ of the sheet, where k is the bending rigidity. 

It is important to note some distinctions in terminology that will be used. The stretching 
(bending) energy given by Us (Ub) denotes the total energy energy contained within an 
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area A in the superlattice. The stretching (bending) energy density is instead denoted by 
lowercase letters us (ub), and can vary with position. The average stretching (bending) 
energy density in the superlattice is Us/A {pis/A). 


B. Elastic Parameters and Superlattice Geometry 

For the numerical results presented here, we use elastic parameters obtained from atom¬ 
istic force field total energy calculations: k = 1.38 eV, and isotropic elastic constants of 
Cxxxx = 26 eV/A^, Cxxyy = 7.5 eV/A^, and Cxyxy = 9.25 eV/A^. These were obtained us¬ 
ing the LAMMPS molecular simulation package^^ and Tersoff potentials'^ with appropriate 
intermixing parameters to describe interfacial interactions. Note that with this choice of 
potentials the T = 0 K lattice constants are ac = 2.46 A and obn = 2.52 A for graphene 
and boron nitride respectively, giving a lattice mismatch of / = (obn — CLc)lcic ~ 2.4%, 
in reasonable agreement with the experimental lattice mismatch of 2.0%^^. We use these 
values for the lattice parameters and the elastic moduli (rather than experimental ones) in 
order to directly compare the continuum results presented here to atomistic simulations. 

In the following analysis, all superlattices are infinitely extended in the plane. To consider 
different wavelength ripples, we denote by [L, H] the number of ‘unit cell’ building blocks 
which are tiled in space to build a graphene/boron nitride supercell into which a single 
wavelength ripple is introduced (along the L direction, parallel to the interface as in Fig. 1). 
Periodic boundary conditions are then imposed in all directions. A single unit cell [L, H] = 
[1,1] is outlined in Fig. 2 which has dimensions 2.46 x 2.46\/3 A^. In our notation, an [L, H] 
superlattice is constructed by repeating this unit cell L times parallel to the interface and 
H times perpendicular to the interface for each subdomain (graphene and boron nitride). 
Thus the area of an [L, H] supercell is given by A = 2Li7(2.46^v^) h?. 


C. Ripple formation in an isolated sheet. 

For simplicity we first consider a thin sheet (say, boron nitride) that has natural (un¬ 
strained) length L(1 -|- /) along the x-direction and natural height H along the y-direction. 
An eigenstrain is introduced, in which the sheet is uniformly compressed in-plane in the 
x-direction: = —/x, = —/. The length of the sheet in the x-direction is now L rather 
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than L{1 + /). The total strain energy in this reference configuration is 

with {CxxxxP/‘^)HL, = 0. (Hereafter, for simplicity we will refer to Cxxxx by 

C so that ~ {Cp/2)HL). 

Ripples of wavelength A will form if the rippled energy U is smaller than . Consider 
(i) the introduction of a ripple w{x) = HA cos(27ra;/A) of wavelength X = L and amplitude 
HA (H is a dimensionless constant to be determined), and (ii) the superposition to the in¬ 
plane displacement field u° of a modulation = (|^) sin (^). The modulation 
helps relieve bending strain concentrated at the peaks and troughs of the ripples as we will 
show. Now the total displacement field is 


Ux = - fx + 
w = HA cos 




The total strain component 


( 3 ) 


^XX 




{dx^x ~\~ dxUx) ~\~ 2 (y^X^^X^ ^x^x^x^x 


/ + AV + /^smy^) 


, , 2 'kx. 
sm" (—) 


( 4 ) 


is largely relieved everywhere when H = v7/^ (neglecting the smaller p contribution). 
This corresponds to the scenario that the total arc length of the rippled system recovers the 
natural, unstrained length of the material: 

■ »> 

where S denotes the elliptic integral of the first kind. One can see from the series expansion 
£’(—4H^7r^) = I (1 -|- ttH^) -|- 0{A^) that the choice H = \pf p satisfies Eq. (5). 

While the choice H = \f]p substantially reduces the stretching energy from the reference 
configuration to a small value lAs ^ ^iP)^ ripple formation comes at the cost of a bending 
energy 

Ub = ~ I [ pV'^w)'^dxdy 
2 Jo Jo 

~ A 

= (4/7r\)y . (6) 
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Thus, whether ripples of a given wavelength A will form depends on a competition between 
the energy = [Cp/2)HX that is offset by relieving the compression and the energy 
U ~ Ub = (4/7r^K)y introduced by the ripple. 

With these expressions, it is possible to roughly estimate a critical wavelength for ripple 
formation in boron nitride due to the lattice mismatch with graphene. By comparing the 
total strain energy densities of the flat vs. rippled systems [Cpp vs. d/vr^/v/A^), the 
critical wavelength is determined to be A ~ 13.2 A. This corresponds to A ~ 5.4 of our unit 
cells, and is found to be in excellent agreement with our atomistic results. In the atomistic 
results, during geometry relaxation we always find that ripples with A < 5 unit cells relax 
back to a fiat configuration, while those with A > 5 unit cells are preserved. 


III. RIPPLE FORMATION IN GRAPHENE - BORON NITRIDE SUPERLAT¬ 
TICES. 

In the superlattices, the lattice mismatch puts the graphene and boron nitride subdomains 
into tension and compression, respectively (see Fig. la). To determine the optimal ripple 
wavelengths, we consider a periodic superlattice that consists of boron nitride {—H < r/ < 0) 
and graphene (0 < y < H) subdomains; the natural lengths of the graphene and boron 
nitride along the interface are L and L(1 + /) respectively (/ now is the lattice mismatch, 
0.024). 

We adopt a reference configuration in which the boron nitride of natural length L(1 + /) 
is kept fiat, but uniformly compressed along the interfacial direction, so that it conforms 
perfectly to the natural graphene length L. This amounts to an eigenstrain = —f and 
corresponding displacement field u° = —fx ioi —H < y < 0 (with all other eigenstrain 
components zero everywhere else). In this reference state, = 0 in graphene, ~ 
{Cp/2)HL in boron nitride, and = 0 everywhere. Thus, the reference total energy 
for the superlattice unit cell of area {2HL) is = {Cp/2)HL and the average energy 
density is U'^^p{2HL) = Cp/4. 

We present two models for rippling in the heterojunctions. In the first, the ‘footprint’ 
of the superlattice along the interfacial direction is fixed at L when ripples form, while in 
the second this constraint is relaxed. For both models, rippling will occur if the energy 
of the rippled configuration U < . Note that the second model is more realistic for 
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freestanding superlattices free to deform in 3D (the first model is primarily used to observe 
trends). Also, for comparison purposes, we note that for superlattices constrained to stay 
flat, ideally the total lattice mismatch / would be evenly accommodated by both materials 
(graphene should be stretched with e^x = //2 and boron nitride squished with exx = ~//2). 
In comparison to the reference configuration, this would reduce the strain energy density to 
Us/{2HL) = Cpj8, which is the optimal distribution of mismatch for flat systems . 


A. Model I - Fixed Footprint. 


Displacements. In the first model, we keep the ‘footprint’ of the rippled superlattice the 
same as in the reference configuration (length L along the interface). To the reference state 
we superpose additional displacement fields, so that the total displacement field is 


w{x,y) 


Uxix,y) 


f ^ 

\1 +expiry)) 

|-/^+ (£)sin 

lo, 



(^), if -H <y<0 
if 0 < y < H 


( 7 ) 


The expression for w corresponds to the introduction of a sinusoidal ripple of wavelength 
X = L and amplitude (Aft^A)/(l + exp{^y)). The parameter /3 sets the decay length of 
the ripple relative to the half-pitch H. If /? is small the ripple amplitude is A^nX in boron 
nitride (y < 0) and decays smoothly to zero across the interface into graphene (y > 0); if P 
is sufficiently small the ripple becomes uniform throughout both subdomains (as /3 —>■ 0, the 
amplitude everywhere becomes AhnX/2). Since boron nitride is in compression but graphene 
is initially strain free in the reference configuration, the former case enables rippling to relieve 
the strain in boron nitride without perturbing graphene. The parameters in the variational 
approach (A(,„, P) are obtained by numerically minimizing the total energy U =Us +Ub- 
Numerical Results. Fig. 3 shows the average formation energy density U / (2HX) computed 
for different rippling wavelengths A and superlattices of varying half-pitch H. From Fig. 3a, 
for a given H, as X increases the average formation energy density initially drops rapidly, 
reaches a minimum, and then begins to increase. The position of the minimum depends on 
the half pitch H. For example when H = 10 unit cells (~ 43 A), the optimal wavelength 
is A ~ 20 unit cells (~ 50 A). When H increases to 400 unit cells (~ 1704 A), the optimal 
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wavelength increases to A ~ 50 unit cells (~ 123 A). Fig. 3b shows the same data, now 
plotted vs. superlattice aspect ratio A/i7, which appears to collapse the data onto a single 
curve beyond the minima (to be discussed later). As shown in the inset of Fig. 3b, in the 
limit of large X/H the strain energy density becomes exactly the same as the optimal flat 
systems. In comparison to our previous atomistic results^^, the trends in Fig. 3 match very 
well to the left of the minima. However, they differ markedly after that: in the atomistic 
models, no observable minimum is found, but instead the average formation energy density 
appears to decay quickly to a flnite, converged value as A increases. 

Contributions to Total Energy U. To understand the trends in Fig. 3, we identify three 
key contributions to the total energy U. The first is the bending energy Ub, evaluated from 
Eq. (2). It is similar to Eq. (6) but now the decay of the ripple across the interface (dw/dy) 
must be included in the Laplacian As before Ub quickly decays as A increases, but 

(as shown later in Eq. (8)) the dw/dy contribution gives rise to an additional term with 
scaling 1/H^. Additionally, there are two dominant (and competing) stretching energy 
contributions: Us ~ + U^g ■ interfacial stretching energy arises from any 

residual strains e^x in the interfacial direction once the ripples have formed. It would be 
zero if a ripple in boron nitride (i) had an amplitude A ~ so that the arc length 

recovers exactly the boron nitride natural length L(l + /), but then (ii) immediately decayed 
across the interface so that the graphene subdomain remains entirely strain free as well. The 
second stretching energy contribution U^g , is the competing term that prevents the above 
scenario. It arises from strains eyy transverse to the interface. The decay of the ripple across 
the interface gives rise to a transverse strain eyy — Fyy = | • The decay constant /5 is 

determined by the competition between (favoring quick decay) and Ug (favoring slow 
decay) near the interface. To assess the relative magnitude of these three contributions, we 
consider the nature of the rippling for small, intermediate, and large \/H. 

Small X/H Behavior. In Eig. 3a, the dotted grey line plots the bending energy density 
Ub approximated simply by Ub/{2HX) where the expression for Ub is directly taken from 
Eq. (6). The 1/A^ decay captures the initial drop in the formation energy density very well 
(particularly for large H), suggesting that for small X/H, most of the strain energy density 
comes from the bending contribution. Our variational results also confirm this: an example 
of a rippled geometry and the associated stretching energy density is shown in Fig. 4a for 
X/H = 0.5, H = 20 unit cells. The geometry shows a completely rippled boron nitride, 


but the ripples quickly decay to zero across the interface into graphene. The corresponding 
stretching energy density distribution us is small and localized to the interfacial region. 
The localized nature can be understood as follows. The ripples in boron nitride completely 
relieve the lattice mismatch (numerically we hnd leaving the interior strain- 

free. The decay constant (3 is such that the ripple decays quickly and the interior of graphene 
is also strain free. By contrast, at the interface where the ripple amplitude is changing, there 
is a residual strain e^x that provides a small but non-zero and a transverse strain eyy 
that provides a small but non-zero . On the whole however, for small X/H the preferred 
geometry is one in which ripples form only in boron nitride to relieve the mismatch strain, 
but then quickly decay into graphene. The residual average strain energy density is mainly 
composed of a bending contribution ub- 

Intermediate X/H Behavior. To explain the increase of the average formation energy 
density in Fig. 3 to the right of the minima, we consider how the stretching energy densities 
vig and vary with aspect ratio X/H. If the ripple nature were to remain the same as 
for small A/iif, the amplitude of the ripple in boron nitride {A^^X — vTA/tt) would increase 
linearly with A, but then still must quickly decay to zero into the graphene. This decay of 
large amplitude wave over a relatively short region would result in an increased transverse 
stretching energy. The outcome of this competition is seen in Fig. 4b: for A/i7 = 5 to 
avoid this energy penalty, the rippling geometry is different. The ripples in boron nitride 
decay somewhat, but persist throughout the graphene subdomain. While this smoothen- 
ing of ripples avoids the large transverse stretching energy penalty, the interior of the two 
subdomains are no longer strain-free (the strain energy is larger now due to residual 
lattice mismatch exx = + Fxx)- Thus, in the intermediate X/H regime (X/H ~ 5), the 

parameters (Ahn,P) are determined as an optimal compromise between U^g and This 

is also evident in Fig. 4b, where the local stretching energy density us is found to be larger, 
and now extends further away from the interface. 

Large X/H Behavior. At large X/H = 50 (Fig. 4c), the ripples are now almost entirely 
uniform throughout both subdomains. This suggests that at large X/H., the ripple geometry 
is dictated by avoiding transverse stretching energy. In fact, in Fig. 4c numerically we find 
that U 5 ~ 0 and vig* ~ Cp/8, mimicking the distribution of the optimal flat super lattices. 
Accordingly, we find numerically that amplitude Ah^ = VJ/tt, and decay /5 —>■ 0 so that the 
lattice mismatch / is shared evenly by a compressive strain of 63 , 3 , = —//2 in boron nitride 
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and a tensile strain oi e^x = f in graphene. (The fixed footprint prevents relaxation to a 
flat configuration, but even so the stretching energy density is the same). This is consistent 
with the inset of Fig. 3b, in which the formation energy density was observed to increase 
back to the value of the flat system. 


B. Model I - Scaling Analysis 


Using our variational expressions for the displacement fields, mathematical expressions 
for the average bending energy density Ub/{2HX) and the average interfacial and transverse 
stretching energy density Wg*/(2H\) and / (2HX) can be obtained. This is useful to 
extract the scaling nature of each contribution. We obtain 

1 


cH 


UbI{2HX) = 


2HX 
1 


A 


niy'^wYdxdy 


-H JO 
\2 


Co + Co 


Uf/{2HX) ={-] {AlC)C, 


Uf^/{2HX) = {At^C)C, + C7{A^, /, /?) 


( 8 ) 

(9) 

( 10 ) 


where each Ci {i = 1, ...5) is a function of 13, and is a non-separable function of (3, f, Ah^- 
The explicit forms of the Ci and T are not provided here - they are straightforward to obtain 
but tedious. Most important are the scaling of the functions with /5, which have leading 
order terms as: Ci ~ /3°, C 2 ~ Us ~ (3^, (74 ~ C^ ~ /d°, J- ~ (3^. The two terms in 
Eq. (10) correspond to contributions from graphene and boron nitride, respectively. Since in 
boron nitride buckling is favorable, is a function that approaches zero when the amplitude 
Abn ~ the decay parameter (3 is large. Relatedly, the function C^ approaches zero 

for sufficiently large (3 (or when Abn ~ 0) since rippling is not favorable in graphene. 

We note several observations from these expressions. The average bending energy density 
has explicit 1/A^ and 1/H^ dependences. By contrast, the average transverse stretching 
energy density has no explicit dependence on L ov H, but only depends on the aspect ratio 
[L/HY. Interestingly, the interfacial stretching term also has no explicit dependence on A 
or H, nor does it depend on aspect ratio X/H. This explains the collapse of the numerical 
results to the right of the minimum in Fig. 3b (where the bending contribution has become 
negligible). For such systems with negligible Ub, the {Abn, (3) fhaf minimize Eq. (9) and 
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(10) can also be numerically determined. We have plotted this line on Fig. 3b as well, and it 
corresponds well to the previous numerical results. The scaling of the two stretching energy 
terms, (i) U*g /{2HX) ^ and (ii) Wg*/{2HX) which comprises of two terms 

and r\j {Aiyn — 2\/J/TTY for graphene and boron nitride respectively, are consistent with 
smoothening of ripples as X/H grows. In the limit X/H oo, Equations (9) and (10) are 
minimized by /5 —>• 0, A^n —t for which U^/{2HX) ~ 0 and /{2HX) (7/^/8, 

recovering again the same strain energy distribution of the optimal flat systems. 


C. Model II - Variable Footprint. 

The first model does not properly capture our atomistic results in which, as the wave¬ 
length A increases, the formation energy drops initially but then appears to reach a steady 
value rather than increase^^. Therefore, we introduce a different variational form for the dis¬ 
placement fields, designed to reduce the transverse stretching energy lA^g in order to maintain 
non-uniform rippling even as the wavelength A grows. For instance, our atomistic simula¬ 
tions show that for A/i7 = 5 as in Fig. 4b, both materials are rippled to large amplitude, 
but the boron nitride subdomain rippling amplitude is somewhat larger to accommodate its 
longer length. 

Therefore, in Model II, we allow the ‘footprint’ of the superlattice to be squeezed from 
its reference value L to some other value L(1 — /*) in the interfacial direction, making it 
advantageous for both the graphene and the boron nitride to form ripples. Consequently, as 
we will show, the difference in the ripple amplitude needed in both materials to relieve lattice 
mismatch is reduced, as is the transverse stretching energy U^g . We find that this delays 
the onset of the increase in the formation energy density with X/H m comparison to Model 
I, and predicts formation energies in very good agreement with our atomistic results. An 
important prediction however is that the formation energy density of the rippled superlattices 
in our atomistic results would eventually also increase as X/H increases. Unfortunately this 
occurs at very large wavelengths that render the supercells computationally prohibitive for 
atomistic simulation - so we are unable to directly verify this prediction. 
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Displacements. In Model II, the variational form of the displacement fields are 



Ac + 


Abn Ac 


1 


A cos 


27rx 

~T~ 


exp (§y) 

(-/ -fs)x+ sin (^), 

-fsX+ (^) sin 


( 11 ) 

( 12 ) 


-H < y < 0 

0 <y < H 

The ripple amplitude now smoothly decays from ^(,nA in boron nitride to AcX in graphene. 
Additionally, the footprint of the super lattice is now X = L{1 — fg), and now there are four 
variational parameters: /« > 0, the amplitudes Ac, Abn, and the decay constant (3. 

With this set of reference displacements, graphene is subjected to an in-plane compressive 
strain of e^x ~ —/s, and boron nitride is subjected to an in-plane compressive strain of 
Cxx ~ —f — fs- The advantage is that amplitudes Ac, Abn can be chosen to optimally relieve 
this in-plane compression, at lower transverse stretching energy cost. The criterion that the 
arc length of each subdomain equal each material’s natural length becomes: 

r-A 



dx — 


2\S{-iAl7,^) 


TT 


— L{1 + f) ^ X{1 + f + fs) 




TT 


(1 + fs) 


(13) 


for boron nitride and graphene respectively. To satisfy this, the amplitudes ~ fs and 
~ / + /«• As fs grows, then Abn and Ac become more similar. The most dramatic 
change to the scaling relations in Eq. (9) for is that now A^^/3^ ^ (Abn — Ac)'^l3'^. 
Thus, in Model II as A/i/ grows it is still possible to maintain large (3 with small transverse 
stretching energy. 

Numerical Results. Figure 5 shows the formation energy densities computed from Model 
II. The behavior is now different from Model I, and both qualitatively and quantitatively 
matches our prior atomistic results^^. The initial sharp drop of the average formation energy 
density with wavelength A still arises from the 1/A^ scaling of the bending energy density 
Ub/{2HL). Now however, the formation energy density drops, remains stable for some time, 
and much more slowly begins to increase. The onset of the increase is slowed by the reduced 
amplitude discrepancy in the two subdomains (Abn ~ V/ + /s/"^ fo Ac ~ \fTslTf), reducing 
the transverse stretching energy. Also, from Fig. 5 we note that the numerically obtained 
data does collapse onto itself at larger Xjli as in model I. The optimal rippling wavelengths 
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again depend on the half-pitch H: when H = 10 nnit cells 43 A), the optimal wavelength 
is A ~ 60 unit cells (~ 150 A). When H increases to 400 unit cells (~ 1700 A), the optimal 
wavelength increases to A ~ 160 unit cells (~ 390 A). 

Analogously in Fig. 6 the ripple geometries for X/H = 0.5,5 and 50, with H = 20 unit 
cells, show first a rippled boron nitride but nearly fiat graphene (as before in Model I), then 
a geometry in which both materials have ripples but a small drop in the ripple amplitude 
across the interface is observable, to finally one where the ripples appear more uniform 
throughout the entire domain (but in fact there remains some amplitude discrepancy, that 
is not perceptible because it is small relative to the total amplitude). In all cases, the 
stretching energy density is small. Although it has begun to grow at A/i7 = 50, it is still 
far from energy density of Model I at this aspect ratio, which had already approached the 
optimal flat value. 

Lastly, we note that even in this model, the aspect ratio \/H cannot increase without 
bound without causing the total energy density to increase. This limit arises from the fact 
that the squishing parameter fg cannot grow arbitrarily large. Since the wave amplitudes 
are approximately given by A^n V f + fs/ TT and Ac ~ then the strain component 

^xx = ^ix + Fxx will become large due to the nonlinear term {dxUxY/2. This itself limits the 
degree of squishing (fg) possible before the stretching energy becomes prohibitively large. 
Numerically, we find that fg never grows beyond 0.14, at which point the lateral stretching 
energy starts to exceed the flat, coherent strain energy. 


IV. CONCLUSIONS 

In summary, a variational continuum model and scaling analysis is developed to describe 
the formation of ripples in planar graphene/boron nitride super lattices. Within periodic 
boundary conditions, superlattices are found to be unstable with respect to formation of 
ripples of appropriate wavelength which efficiently accommodate lattice mismatch strain. A 
competition between the bending energy (favoring large wavelength ripples) and transverse 
stretching energy (favoring small wavelength ripples) gives rise to a set of optimal rippling 
wavelengths that minimize the total energy of the deformed superlattice. These optimal 
wavelengths depend on the superlattice pitch, but vary from 15-40 nm for the geometries 
considered here. In principal, the predictions presented may be validated in experiment for 
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freestanding superlattices free to deform out of plane. 
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FIG. 1. (a) A flat planar superlattice of composed of regular, alternating domains of graphene 
(black) and boron nitride (red). Due to the mismatch in the lattice constants (2.46 A and 2.52 
A, respectively), when atomic registry is maintained across the interfaces the flat superlattice is 
composed of alternating regions of tension and compression in the direction parallel to the interface, 
(b) Much of the lattice mismatch strain can be accomodated by the formation of a ripple network, 
in which the ripples of larger amplitude in boron nitride decay to ripples of smaller amplitude in 
graphene. 



FIG. 2. A single [L, i7] = [1,1] unit cell is outlined in red on the honeycomb lattice. The dimensions 
are 2.46 A, 2.46\/3 A for L^H respectively. 
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3.6 5.4 

aspect ratio X/H 


FIG. 3. (a) Average formation energy density for Model I rippled superlattices as a function of 
ripple wavelength A for different half-pitch H. For all i7, initially the formation energy density 
drops rapidly due to the ^ 1/A^ nature of the bending energy. As the wavelength increases further, 
the transverse stretching energy due to decay of ripples across the interface becomes dominant, 
inducing a transition to uniform ripples and an increase in formation energy density, (b) Plotting 
the formation energy density for Model I rippled superlattices as a function of aspect ratio X/H 
collapses the curves together in the large X/H regime. (Inset: as X/H further grows, the formation 
energy density increases to the value of the coherent, flat superlattices). 
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a) X/H = 0,5 (b) X/H = 5 (c) X/H = 50 



0+0 3.6 


FIG. 4. Rippling structure and stretching energy density us of Model I superlattices with half-pitch 
i7=20 unit cells, for aspect ratios X/H = 0.5, 5, and 50. For small X/H the ripples are present 
in boron nitride and decay into graphene. At intermediate A/i7, rippling becomes somewhat 
preserved across the graphene subdomain. The stretching energy density us is no longer localized 
to the interface. For the largest X/H the ripple amplitude is almost uniform across the two 
materials, and is such that the mismatch strain is equally shared by them. In this limit, the boron 
nitride is in compression and the graphene in tension, very similar to the flat superlattices, and 
the formation energy density is the same as in the flat case. 
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FIG. 5. (a) Average formation energy density for Model II rippled superlattices as a function of 
ripple wavelength A = L(1 — fg) for different half-pitch H. For all i7, (as in Model I) initially the 
formation energy density drops rapidly due to the ^ 1/A^ nature of the bending energy. However, 
since both materials are now compressed along the interfacial direction, ripple formation is possible 
in both. The reduced amplitude discrepancy between and Ac effectively delays the onset of 
the increase in the formation energy density with increasing X/H. (b) The same data plotted as 
a function of aspect ratio X/H] again the results collapse to the same curve for sufficiently large 
aspect ratio. 
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FIG. 6. Rippling structure and stretching energy density us of Model II superlattices with half¬ 
pitch H=20 unit cells, for aspect ratios X/H = 0.5, 5, and 50. The stretching energy density us 
is very small for all geometries shown (and largely composed of numerical noise - we maintain the 
same scale bar as in Fig. 4, to facilitate comparison). 
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